% 这个脚本演示了二维柱面波的叠加，使用了解析解与远场近似
% 在波程差为波长整数倍的地方可以看到相长干涉等
% CC-BY 4

clc
clear

[x y]=meshgrid(-5:0.1:5,-5:0.1:5);
[nx ny] = size(x);
z = zeros(nx,ny);

v = 1;

source(1).x=-2;
source(1).y=0;
source(1).T=1;
source(1).w = 2*pi/source(1).T;
source(1).k = source(1).w/v;

source(2).x=2;
source(2).y=0;
source(2).T=1;
source(2).w = 2*pi/source(2).T;
source(2).k = source(2).w/v;

t = 0;
picind = 1;

figure();

for t = 0:0.1:20

  for i = 1:nx
    for j = 1:ny
      u = 0;
      dist1 = sqrt((x(i,j)-source(1).x)^2 + (y(i,j)-source(1).y)^2);
      dist2 = sqrt((x(i,j)-source(2).x)^2 + (y(i,j)-source(2).y)^2);

      u += (1/dist1)*sin(source(1).w*t-source(1).k*dist1);
      u += -(1/dist2)*sin(source(2).w*t-source(2).k*dist2);

      if isnan(u) || isinf(u) || abs(u) > 2
        if u > 0
          u = 2;
        else
          u = -2;
        endif
      endif

      z(i,j) = u;
    endfor
  end

  clf;
  hold on
  axis equal
  axis([-5 5 -5 5 -5 5])
  f = surf(x,y,z);
  caxis([-1 1])
  set(f,'edgecolor','none');

  drawnow
  pause(0.01);
##  print(['wave_2D' num2str(picind) '.png'],'-dpng')
  picind++;
end

